# preamble #### 
rm(list=ls())
## set your working directory to the root folder for this replication file here
setwd('~/Dropbox/CATF EU Survey/replication')
figures<-"figures/"
tables<-"tables/"

library(estimatr)
library(randomizr)
library(DeclareDesign)
library(haven) ## to read sav files
library(WeightIt)
library(texreg)
library(car)
library(emmeans)
library(ggpubr)
library(xtable)
library(RColorBrewer)
library(scales)
library(magrittr)
library(tidyverse)
set.seed(1234)

# load data #### 
load("BergquistMahdavi_data.Rda")
## create unique id for each person 
d<-d%>%
  mutate(id=paste(country,RecordNo,sep="-"))

# recode demographic variables #### 
## relevel ideology variable
d$ideo<-relevel(factor(d$ideo),ref="other")




## scale education and ideology variables within each country 
## voceduc_neu is a 5-point code for germany: 
# <1> No degree
# <2> Still in training
# <3> Still studying
# <4> Apprenticeship or comparable qualification (secondary or below)
# <5> University or technical college degree (post-secondary)
# <777> not specified

## education_level_it in Italy 
# <1> No diploma 
# <2> Elementary school certificate 
# <3> Middle school diploma 
# <4> High school diploma 
# <5> Apprenticeship diploma 
# <6> Technical Institute (secondary school or below)
# <7> Bachelor's degree (technical degree or college)
# <8> Single cycle master's degree (technical degree or college)
# <9> Old system degree
# <10> First level master's degree
# <11> Specialist degree/ master's degree
# <12> Second level doctorate / master's degree
# <13> Other
# <14> I prefer not to answer

## education is a 7-category scale in France
# <1> No degree--1
# <2> Study certificate--below high school 
# <3> College Certificate--below high school (get this at age 15)
# <4> Short technical education diploma (CAP, BEP or equivalent) (less than high school)
# <5> Baccalaureate, professional certificate or equivalent (this is equivalent to high school)
# <6> Undergraduate or technical diploma (Bac +2, BTS, DEUG, License) (post-secondary technical degree) 
# <7> Second/third university degree (Master, Grandes Ecoles, Doctorate) (academic post-secondary degree)
# <8> Other
# <9> I don't know / I prefer not to answer

## education_level2_pl is a 7-category scale in Poland
# <1> Elementary/Middle School
# <2> Secondary/Post-secondary
# <3> Professional
# <4> Higher
# <5> Doctorate/Higher degree
# <6> I prefer not to answer

d<-d%>%
  mutate(education_recode=
           case_when(country=="Germany"&voceduc_neu==5~"Post-secondary degree",
                     country=="Germany"&voceduc_neu<5~"Secondary degree or less",
                     country=="Italy"&education_level_it<7~"Secondary degree or less",
                     country=="Italy"&education_level_it>=7&education_level_it<=12~"Post-secondary degree",
                     country=="France"&education<=5~"Secondary degree or less",
                     country=="France"&education>5&education<8~"Post-secondary degree",
                     country=="Poland"&education_level2_pl<=2~"Secondary degree or less",
                     country=="Poland"&education_level2_pl>2&education_level2_pl<6~"Post-secondary degree"
           ))
table(d$education_recode,d$country)

## recode income 
# germany: hinc (monthly)
# france: income (annual)
# italy: profile_gross_household_EU (annual)
# poland: income_PL () (monthly)

## recode to midpoint of ranges and then standardize within countries
d<-d%>%
  mutate(income_recode=case_when(income_PL==1~2000,
                                 income_PL==2~2500,
                                 income_PL==3~3500,
                                 income_PL==4~5000,
                                 income_PL==5~7000,
                                 income_PL==6~8000,
                                 hinc==1~500,
                                 hinc==2~750,
                                 hinc==3~1250,
                                 hinc==4~1750,
                                 hinc==5~2250,
                                 hinc==6~2750,
                                 hinc==7~3250,
                                 hinc==8~3750,
                                 hinc==9~4250,
                                 hinc==10~4750,
                                 hinc==11~7500,
                                 hinc==12~10000,
                                 income==1~15000,
                                 income==2~17500,
                                 income==3~25000,
                                 income==4~35000,
                                 income==5~45000,
                                 income==6~55000,
                                 income==7~65000,
                                 income==8~75000,
                                 income==9~85000,
                                 income==10~95000,
                                 income==11~112500,
                                 income==12~137500,
                                 income==13~175000,
                                 income==14~200000,
                                 profile_gross_household_EU==1~5000,
                                 profile_gross_household_EU==2~7500,
                                 profile_gross_household_EU==3~12500,
                                 profile_gross_household_EU==4~17500,
                                 profile_gross_household_EU==5~22500,
                                 profile_gross_household_EU==6~27500,
                                 profile_gross_household_EU==7~32500,
                                 profile_gross_household_EU==8~37500,
                                 profile_gross_household_EU==9~42500,
                                 profile_gross_household_EU==10~47500,
                                 profile_gross_household_EU==11~55000,
                                 profile_gross_household_EU==12~65000,
                                 profile_gross_household_EU==13~80000,
                                 profile_gross_household_EU==14~125000,
                                 profile_gross_household_EU==15~325000,
                                 profile_gross_household_EU==16~500000
  ))
## standardize income within countries
standardize<-function(x)(x-mean(x,na.rm=TRUE))/sd(x,na.rm=TRUE)
d<-d%>%
  group_by(country)%>%
  mutate(incomestd=standardize(income_recode))%>%
  ungroup()

# Adding an age group category for all 
d <- d %>% 
  mutate(
    age_group_all = case_when(
      age < 25 ~ "18-24",
      age %in% 25:34 ~ "25-34",
      age %in% 35:44 ~ "35-44",
      age %in% 45:54 ~ "45-54",
      age %in% 55:64 ~ "55-64",
      age > 64 ~ "65+"
    )
  )


# Table S1: demographic crosstabs of our DV #### 
# note that we formatted this table manually in LaTex
# Summaries for each group
desc.all <- d %>% 
  group_by(country) %>%
  summarize(
    group = "all",
    pctsupport = mean(dv.methanelaw.binary, na.rm = T)
  )

desc.age <- d %>%
  group_by(country,age_group_all) %>%
  summarize(
    pctsupport = mean(dv.methanelaw.binary, na.rm = T)
  )

desc.ideo <- d %>%
  group_by(country,ideo) %>%
  summarize(
    pctsupport = mean(dv.methanelaw.binary, na.rm = T)
  )

desc.educ <- d %>%
  filter(!is.na(education_recode)) %>%
  group_by(country,education_recode) %>%
  summarize(
    pctsupport = mean(dv.methanelaw.binary, na.rm = T)
  ) 

desc.gend <- d %>%
  group_by(country,gender) %>%
  summarize(
    pctsupport = mean(dv.methanelaw.binary, na.rm = T)
  ) 

names(desc.age) <- names(desc.ideo) <- names(desc.educ) <- names(desc.gend) <- c("country","group","pctsupport")

desc.table <- data.frame(rbind(desc.all,desc.age,desc.ideo,desc.educ,desc.gend)) %>%
  spread(country, pctsupport) %>%
  arrange(group)

print(xtable(desc.table),include.rownames=FALSE)

# Tables S2, S3: are the country-level responses to our DV significantly different from each other? ####
reg1<-lm_robust(dv.methanelaw.binary~country,data=d)
summary(reg1)

texreg(reg1,include.ci = FALSE,#custom.model.names=c("LPM Model","Linear model with ordinal DV"),
       caption="Differences in policy support between countries",
       file=paste0(tables,"countrydiffs.tex"),
       label="tab:countrydiffs",float.pos="!htpb")

lh1<-tidy(linearHypothesis(reg1,"countryGermany=countryItaly"))
lh2<-tidy(linearHypothesis(reg1,"countryGermany=countryPoland"))
lh3<-tidy(linearHypothesis(reg1,"countryItaly=countryPoland"))
lh<-do.call(rbind,list(lh1,lh2,lh3))%>%
  mutate(term=paste0(term," = 0"))%>%
  rename(`Hypothesis test`=term)%>%
  select(`Hypothesis test`,statistic,p.value)%>%
  rename(`Chi-squared statistic`=statistic)
print.xtable(xtable(lh,caption="Wald tests of equality of coefficients for country effects on responses to the dependent variable used in our experiments",
                    label="tab:countrydiffs.lh"),include.rownames =FALSE,file=paste0(tables,"countrydiffs_lhtests.tex"))

# table s4: differences by ideology within country #### 
d$ideo<-relevel(d$ideo,ref="right")
reg2<-lm_robust(dv.methanelaw.binary~ideo,subset(d,country=="France"))
reg3<-lm_robust(dv.methanelaw.binary~ideo,subset(d,country=="Germany"))
reg4<-lm_robust(dv.methanelaw.binary~ideo,subset(d,country=="Italy"))
reg5<-lm_robust(dv.methanelaw.binary~ideo,subset(d,country=="Poland"))
texreg(list(reg2,reg3,reg4,reg5),include.ci=FALSE,
       custom.model.names=c("France","Germany","Italy","Poland"),
       caption="Differences in policy support between ideology groups within countries. The table shows the effect of indicators for each ideological grouping (proxied using vote choice in the most recent national election) on support for methane regulations. The indicator for belonging to the ideological right is the omitted category.",
       file=paste0(tables,"ideodiffs.tex"),
       label="tab:ideodiffs",float.pos="!htpb")

# sample size across countries #### 
table(d$country)

# table S5: covariate balance ####
## make gender, ideo, education numeric
d.bal <- d %>%
  mutate(
    gender.n = case_when(gender=="Male"~0,
                         gender=="Female"~1),
    ideo.n = case_when(ideo=="left"~-1,
                       ideo=="other"~0,
                       ideo=="right"~1),
    educ.n = case_when(education_recode=="Secondary degree or less"~0,
                       education_recode=="Post-secondary degree"~1)    
  )

## ANOVAs - Cost
age.aov.cost <- aov(age ~ z.cost, data = d.bal)
gen.aov.cost <- aov(gender.n ~ z.cost, data = d.bal)
ide.aov.cost <- aov(ideo.n ~ z.cost, data = d.bal)
inc.aov.cost <- aov(incomestd ~ z.cost, data = d.bal)
edu.aov.cost <- aov(educ.n ~ z.cost, data = d.bal)

## ANOVAs - Frame
age.aov.frame <- aov(age ~ z.frame, data = d.bal)
gen.aov.frame <- aov(gender.n ~ z.frame, data = d.bal)
ide.aov.frame <- aov(ideo.n ~ z.frame, data = d.bal)
inc.aov.frame <- aov(incomestd ~ z.frame, data = d.bal)
edu.aov.frame <- aov(educ.n ~ z.frame, data = d.bal)


## Tukey's multiple comparisons test - Cost
age.thsd.cost <- tidy(TukeyHSD(age.aov.cost))
gen.thsd.cost <- tidy(TukeyHSD(gen.aov.cost))
ide.thsd.cost <- tidy(TukeyHSD(ide.aov.cost))
inc.thsd.cost <- tidy(TukeyHSD(inc.aov.cost))
edu.thsd.cost <- tidy(TukeyHSD(edu.aov.cost))

## Tukey's multiple comparisons test - Frame
age.thsd.frame <- tidy(TukeyHSD(age.aov.frame))
gen.thsd.frame <- tidy(TukeyHSD(gen.aov.frame))
ide.thsd.frame <- tidy(TukeyHSD(ide.aov.frame))
inc.thsd.frame <- tidy(TukeyHSD(inc.aov.frame))
edu.thsd.frame <- tidy(TukeyHSD(edu.aov.frame))

rm(age.aov.cost, gen.aov.cost, ide.aov.cost, inc.aov.cost, edu.aov.cost,
   age.aov.frame, gen.aov.frame, ide.aov.frame, inc.aov.frame, edu.aov.frame)

# As data frame
age.thsd.cost$variable <- "age"
gen.thsd.cost$variable <- "gender"
ide.thsd.cost$variable <- "ideology"
inc.thsd.cost$variable <- "income"
edu.thsd.cost$variable <- "education"

age.thsd.frame$variable <- "age"
gen.thsd.frame$variable <- "gender"
ide.thsd.frame$variable <- "ideology"
inc.thsd.frame$variable <- "income"
edu.thsd.frame$variable <- "education"

thsd <- rbind(age.thsd.cost, gen.thsd.cost, ide.thsd.cost, inc.thsd.cost, edu.thsd.cost,
              age.thsd.frame, gen.thsd.frame, ide.thsd.frame, inc.thsd.frame, edu.thsd.frame)

thsd <- thsd %>% select(term, contrast, variable, estimate, adj.p.value) %>% arrange(term,contrast,variable)
bal.xt <- xtable(thsd, 
                 digits = 3, 
                 caption = "Covariate balance table using ANOVA with p-values adjusted using Tukey's multiple comparison test.",
                 label = "tab:balance")

names(bal.xt) <- c("Treatment","Condition comparison","Variable","Difference in means","Tukey-adjusted p-value")
print(bal.xt, include.rownames=FALSE, file=paste0(tables,"results_balance.tex"))

# descriptive figures from survey #### 
## Fig 1: associations with methane #### 
dphrases<-d%>%
  select(id,country,contains("phrasesassociate"),NG_Methane)
dphrases<-dphrases%>%
  pivot_longer(contains("phrase"),values_to='val',names_to='phrase')
dphrases<-dphrases%>%
  mutate(phrase=parse_number(phrase),
         phrase=case_when(phrase==1~"Cows and other animals",
                          phrase==2~"Heating and cooking",
                          phrase==3~"Greenhouse gas",
                          phrase==4~"Clean burning fuel",
                          phrase==5~"Air pollution",
                          phrase==6~"Fracking",
                          phrase==7~"Permafrost melting",
                          phrase==8~"Fossil fuel",
                          phrase==9~"Transition fuel",
                          phrase==10~"Energy security"))
dphrases<-dphrases%>%filter(val!=0)
prop.table(table(dphrases$phrase))*100
## percent of people who chose each phrase
dphrases.sum<-dphrases%>%
  group_by(phrase,NG_Methane)%>%
  summarise(val=sum(val))%>%
  mutate(prop=val*100/length(unique(dphrases$id)))%>%
  mutate(NG_Methane=ifelse(NG_Methane==1,"Associations with\nNatural gas","Associations with\nMethane"))#%>%

dphrases.sum[dphrases.sum$NG_Methane=="Associations with\nMethane",]
phraseorder<-dphrases.sum[dphrases.sum$NG_Methane=="Associations with\nMethane",]%>%arrange(prop)
phraseorder<-phraseorder$phrase
dphrases.sum$phrase<-factor(dphrases.sum$phrase,levels=phraseorder)

dphrases.sum.country<-dphrases%>%
  group_by(phrase,country,NG_Methane)%>%
  summarise(val=sum(val))%>%
  mutate(prop=val*100/length(unique(dphrases$id)))%>%
  mutate(NG_Methane=ifelse(NG_Methane==1,"Natural gas","Methane"))
dphrases.sum.country$phrase<-factor(dphrases.sum.country$phrase,levels=phraseorder)
dphrasefig.country<-dphrases.sum.country%>%
  ggplot(aes(x=prop,y=phrase,fill=NG_Methane))+
  theme_bw()+
  geom_bar(stat="identity",position="dodge")+
  scale_fill_brewer(palette="Dark2",guide=guide_legend(reverse=TRUE))+
  labs(x="Percent of subjects who chose each phrase",
       y="Phrases associated with\n`methane' and `natural gas'",fill="")+
  theme(legend.position="bottom")+
  scale_x_continuous(limits=c(0,10),breaks=c(0,5,10))+
  facet_wrap(~country,ncol=4,nrow=1)
dphrasefig.country
ggsave(file=paste0(figures,"descriptive_phraseassociations_bycountry.pdf"),width=6.5,height=3.5,dphrasefig.country)

## questions about methane composition in the atmosphere #### 
## Fig. S1: knowledge about methane as a component of natural gas #### 
gascomptab<-data.frame(table(d$gascomposition,d$country))
gascompn<-gascomptab%>%
  group_by(Var2)%>%
  summarise(n=sum(Freq))
gascomptab<-gascomptab%>%left_join(gascompn,by="Var2")%>%
  mutate(prop=Freq*100/n)%>%
  select(Var2,Var1,prop)%>%
  mutate(Var1=factor(Var1,levels=c("don't know","Carbon dioxide","Hydrogen","Oxygen","Methane")))
fig1<-gascomptab%>%
  ggplot(aes(x=prop,y=Var2,fill=Var1))+
  theme_bw()+
  geom_bar(stat="identity",position="stack")+
  # scale_fill_brewer(palette="Dark2",guide=guide_legend(reverse=TRUE))+
  scale_fill_manual(values=brewer.pal(8,"Dark2")[c(8,1:4)],guide=guide_legend(reverse=TRUE))+
  labs(x="Percent of respondents",y="",fill="Which composes\nthe largest percentage\nof natural gas?")+
  theme(legend.position="left")
fig1

gascompfig<-
  d%>%
  ggplot(aes(x=gasmethane))+
  geom_density()+
  geom_area(aes(x=stage(gasmethane,after_stat=oob_censor(x,c(69.9,94.9)))),
            stat="density",fill=brewer.pal(8,"Dark2")[1],alpha=.5)+
  labs(x="Methane as percent of natural gas composition",y="Density of responses")+
  theme_classic()+
  theme(axis.text.y=element_blank(),
        axis.ticks.y=element_blank())

## number of responses in the "correct" range 
nrow(d[d$gasmethane>69.9&d$gasmethane<95,])/nrow(d)

ggsave(paste0(figures,"descriptive_methanegascomp.pdf"),width=6.5,height=6.5,
       ggarrange(fig1+theme(legend.position="left"),gascompfig,nrow=2,ncol=1))

## Fig. 2: knowledge, perceptions, concern about methane #### 
d<-d%>%
  mutate(gascomposition.binary=ifelse(gascomposition=="Methane",1,0))

### responses to change in methane in atmosphere question #### 
methanechangetab<-data.frame(table(d$methanechange,d$country))
methanechangen<-methanechangetab%>%
  group_by(Var2)%>%
  summarise(n=sum(Freq))
methanechangetab<-methanechangetab%>%left_join(methanechangen,by="Var2")%>%
  mutate(prop=Freq*100/n)%>%
  select(Var2,Var1,prop)%>%
  mutate(Var1=factor(Var1,levels=c("stayed same","decreased","increased")))
fig2<-methanechangetab%>%
  ggplot(aes(x=prop,y=Var2,fill=Var1))+
  theme_bw()+
  geom_bar(stat="identity",position="stack")+
  # scale_fill_brewer(palette="Dark2",guide=guide_legend(reverse=TRUE))+
  scale_fill_manual(values=brewer.pal(8,"Dark2")[2:4],guide=guide_legend(reverse=TRUE))+
  labs(x="Percent of respondents",y="",fill="Has methane in the\natmosphere increased,\ndecreased, or stayed\nthe same?")+
  theme(legend.position="left")
fig2

### question about methane as a problem ####
probtab<-data.frame(table(d$country,d$methaneproblem))
probn<-probtab%>%
  group_by(Var1)%>%
  summarise(n=sum(Freq))%>%
  ungroup()
probtab<-left_join(probtab,probn,by="Var1")%>%
  mutate(prop=Freq*100/n)%>%
  select(Var1,Var2,prop)
probfig<-probtab%>%
  ggplot(aes(x=prop,y=Var1,fill=Var2))+
  theme_bw()+
  geom_bar(stat="identity",position="stack")+
  # scale_fill_brewer(palette="Dark2",guide=guide_legend(reverse=TRUE))+
  scale_fill_manual(values=brewer.pal(8,"Dark2")[c(8,2:4)],guide=guide_legend(reverse=TRUE))+
  labs(x="Percent of respondents",y="",fill="How much of a problem\nis methane\nfor the climate?")+
  theme(legend.position="left")
probfig

# save figure as two panels
ggsave(paste0(figures,"descriptive_methaneknowledge_beliefs.pdf"),width=6.5,height=6.5,
       ggarrange(fig2+labs(x=""),probfig,
       ncol=1,nrow=2))




# Fig 3: policy questions that precede our experiment ####
## collapse methane fee questions (cost info doesn't make much of a difference, so note in paper that this marginalizes across random assignment of cost and benefit info)
fee<-d%>%
  select(importregulations_2a,importregulations_3a,importregulations_4a)%>%
  pivot_longer(everything(),names_to="var",values_to="val")%>%
  filter(!is.na(val))
feesum<-data.frame(prop.table(table(fee$val))*100)%>%mutate(var="Impose fee on\nimported oil and gas")
importregs<-data.frame(prop.table(table(d$importregulations_1a))*100)%>%mutate(var="Regulate imported gas")
emissionregs<-data.frame(prop.table(table(d$methaneemissiona))*100)%>%mutate(var="Regulate methane leaks")
regcompanies<-data.frame(prop.table(table(d$financialincentive))*100)%>%mutate(var="Regulate companies")
techupdate<-data.frame(prop.table(table(d$techupdatea))*100)%>%mutate(var="Require equipment upgrades")
regfig<-do.call(rbind,list(feesum,importregs,emissionregs,techupdate))%>%
  mutate(response=case_when(Var1==0~"Strongly oppose",
                            Var1==1~"Oppose",
                            Var1==2~"Support",
                            Var1==3~"Strongly support"),
         response=factor(response,levels=c("Strongly support","Support","Oppose","Strongly oppose")))%>%
  ggplot(aes(x=Freq,y=var,fill=response))+
  theme_bw()+
  geom_bar(stat="identity",position="stack")+
  scale_fill_brewer(palette="Dark2",guide=guide_legend(reverse=TRUE))+
  labs(x="Percent of respondents",y="Policy proposal",fill="Response")
regfig
ggsave(file=paste0(figures,"descriptive_policysupport.pdf"),width=6.5,height=3.5,regfig)

# Experimental analysis  ####
## cost sensitivity regressions #### 
# set up IPW weights 
sum(table(d$z.cost))
w <- weightit(z.cost~block,data=d,estimand="ATE",method="ps")
d$w.cost<-w$weights
d$z.cost<-relevel(factor(d$z.cost),ref="no cost info")
h1mod<-lm_robust(dv.methanelaw.binary~z.cost+country+ideo+gender+age+education_recode+income_recode,weights=w.cost,data=d) 

## apply Benjamini-Hochberg correction
bh.h1mod <- round(p.adjust(h1mod$p.value[grepl("z",names(h1mod$p.value))],"BH"),2)

## analysis with ordinal outcome as continuous variable 
h1mod.cont<-lm_robust(dv.methanelaw~z.cost+country+ideo+gender+age+education_recode+income_recode,weights=w.cost,data=d)
bh.h1mod.cont <- round(p.adjust(h1mod.cont$p.value[grepl("z",names(h1mod$p.value))],"BH"),2)

## Table S6: cost results
texreg(list(h1mod,h1mod.cont),include.ci = FALSE,custom.model.names=c("LPM Model","Linear model with ordinal DV"),
       caption=paste("Effect of cost information on policy support. Benjamini-Hochberg adjusted p-values for cost coefficients:", bh.h1mod[1],  bh.h1mod[2], "(LPM)",bh.h1mod.cont[1],bh.h1mod.cont[2],"(Continuous DV)", sep = " "),
       file=paste0(tables,"costresults.tex"),
       label="tab:costresults",float.pos="!htpb")

## benefit frames regressions #### 
w<-weightit(z.frame~block,estimand="ATE",method="ps",data=d)
d$w.frame<-w$weights
d$z.frame<-relevel(factor(d$z.frame),ref="no frame")
h2mod<-lm_robust(dv.methanelaw.binary~z.frame+country+ideo+gender+age+education_recode+income_recode,weights=w.frame,data=d) 
## extract Benjamini-Hochberg corrected p values
bh.h2mod <- round(p.adjust(h2mod$p.value[grepl("z.frame",names(h2mod$p.value))],"BH"),2)

## continuous DV
h2mod.cont<-lm_robust(dv.methanelaw~z.frame+country+ideo+gender+age+education_recode+income_recode,weights=w.frame,data=d) 
## extract Benjamini-Hochberg corrected p values
bh.h2mod.cont <- round(p.adjust(h2mod.cont$p.value[grepl("z.frame",names(h2mod$p.value))],"BH"),2)

## table S7: framing results
texreg(list(h2mod,h2mod.cont),include.ci=FALSE,custom.model.names = c("LPM Model","Linear model with ordinal DV"),
       caption=paste("Effects of benefits framing on policy support. Benjamini-Hochberg adjusted p-values for benefit coefficients:", bh.h2mod[1], bh.h2mod[2], bh.h2mod[3], "(LPM model)",
                     bh.h2mod.cont[1], bh.h2mod.cont[2], bh.h2mod.cont[3], "(Linear model with ordinal DV)",sep = " "),
       file=paste0(tables,"frameresults.tex"),
       label="tab:frameresults",float.pos="!htpb")

## cost x benefit regressions #### 
## code interacted treatment groups
table(d$z.cost)
table(d$z.frame)
d<-d%>%
  mutate(z=case_when(z.cost=="no cost info"&z.frame=="no frame"~"control",
                     z.cost=="high"&z.frame=="climate"~"increase.climate",
                     z.cost=="low"&z.frame=="climate"~"decrease.climate",
                     z.cost=="high"&z.frame=="energysecurity"~"increase.security",
                     z.cost=="low"&z.frame=="energysecurity"~"decrease.security",
                     z.cost=="high"&z.frame=="health"~"increase.health",
                     z.cost=="low"&z.frame=="health"~"decrease.health"))

w<- weightit(z~block,data=d,estimand="ATE",method="ps")
d$w.full<-w$weights

h3amod<-lm_robust(dv.methanelaw.binary~z+country+ideo+gender+age+education_recode+income_recode,
                  data=d,weights=w.full)

h3amod.cont<-lm_robust(dv.methanelaw~z+country+ideo+gender+age+education_recode+income_recode,
                       data=d,weights=w.full)

## table S8: cost sensitivity of benefits framing
texreg(list(h3amod,h3amod.cont),
       custom.model.names=c("LPM Model","Linear model with ordinal DV"),
       custom.coef.names = c("(Intercept)",
                             "Cost.decrease:Benefit.climate",
                             "Cost.decrease:Benefit.health",
                             "Cost.decrease:Benefit.security",
                             "Cost.increase:Benefit.climate",
                             "Cost.increase:Benefit.health",
                             "Cost.increase:Benefit.security",
                             "Germany","Italy","Poland","Left","Right","Male",
                             "Age","Secondary degree or lower","Income"),
       include.ci=FALSE,caption="Cost sensitivity of benefit framing",
       file=paste0(tables,"results_interactions.tex"),
       label="tab:interactions",float.pos="!htpb")

## wald test for cost-sensitivity of response to benefit frames #### 
waldclimate<-as.data.frame(linearHypothesis(h3amod,"zdecrease.climate=zincrease.climate"))[2,]%>%
  mutate(`Benefit frame`="Climate")%>%
  select(`Benefit frame`,`Chisq`,`Pr(>Chisq)`)

waldhealth<-as.data.frame(linearHypothesis(h3amod,"zdecrease.health=zincrease.health"))[2,]%>%
  mutate(`Benefit frame`="Health")%>%
  select(`Benefit frame`,`Chisq`,`Pr(>Chisq)`) ## this one is significant, but in the wrong direction. (ie, effect in the cost-increasing condition is positive, but null. effect in cost-decreasing condition is negative, and significant)

waldsecurity<-as.data.frame(linearHypothesis(h3amod,"zdecrease.security=zincrease.security"))[2,]%>%
  mutate(`Benefit frame`="Security")%>%
  select(`Benefit frame`,`Chisq`,`Pr(>Chisq)`)
wald.tests<-do.call(rbind,list(waldclimate,waldhealth,waldsecurity))

## table S9: interactions between cost and benefit frames #### 
print.xtable(xtable(wald.tests,caption="Wald tests for interaction effects between cost and benefit frames",
                    label="tab:waldtests"),include.rownames = FALSE,file=paste0(tables,"results_waldtests.tex"))

linearHypothesis(h3amod.cont,"zdecrease.climate=zincrease.climate")
linearHypothesis(h3amod.cont,"zdecrease.health=zincrease.health")
linearHypothesis(h3amod.cont,"zdecrease.security=zincrease.security")

## 2-panel plot with marginalized and fully interacted effects #### 
frameresults<-tidy(h2mod)
costresults<-tidy(h1mod)
intresults<-tidy(h3amod)

frameresults.cont<-tidy(h2mod.cont)
costresults.cont<-tidy(h1mod.cont)
intresults.cont<-tidy(h3amod.cont)

results<-rbind(frameresults,costresults)

results.cont<-rbind(frameresults.cont,costresults.cont)

mainplot<-results%>%
  filter(grepl("z.",term))%>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error,
         ci.high95=estimate+1.96*std.error,
         ci.low95=estimate-1.96*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_linerange(aes(xmin=ci.low95,xmax=ci.high95),lwd=.15)+
  geom_vline(xintercept=0,lty="dotted")+
  theme_classic()+
  scale_y_discrete(
    labels = c("z.costlow" = "Lower costs",
               "z.costhigh" = "Higher costs",
               "z.framehealth" = "Health",
               "z.frameclimate" = "Climate change",
               "z.frameenergysecurity" = "Energy security"))+
  labs(x="",y="")+
  coord_cartesian(xlim = c(-0.1,0.05)) +
  theme(axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=11))
mainplot

mainplot.cont<-results.cont%>%
  filter(grepl("z.",term))%>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  theme_classic()+
  scale_y_discrete(
    labels = c("z.costlow" = "Lower costs",
               "z.costhigh" = "Higher costs",
               "z.framehealth" = "Health",
               "z.frameclimate" = "Climate change",
               "z.frameenergysecurity" = "Energy security"))+
  labs(x="",y="")+
  coord_cartesian(xlim = c(-0.2,0.1)) +
  theme(axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=11))
mainplot.cont

intplot<-intresults%>%
  filter(grepl("z.",term))%>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error,
         ci.high95=estimate+1.96*std.error,
         ci.low95=estimate-1.96*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_linerange(aes(xmin=ci.low95,xmax=ci.high95),lwd=.15)+
  geom_vline(xintercept=0,lty="dotted")+
  scale_y_discrete(limits=c("zincrease.climate","zdecrease.climate","zincrease.security","zdecrease.security","zincrease.health","zdecrease.health"),
                   labels = c("zincrease.climate"= "Higher cost,\nclimate change","zdecrease.climate"="Lower cost,\nclimate change", 
                              "zincrease.security"="Higher cost,\nenergy security","zdecrease.security"="Lower cost,\nenergy security",
                              "zincrease.health"="Higher cost,\nhealth","zdecrease.health"="Lower cost,\nhealth")
  ) +
  labs(x="",y="")+
  coord_cartesian(xlim = c(-0.1,0.05)) +
  theme_classic()+
  theme(axis.text.x = element_text(size=11), 
        axis.text.y = element_text(size=11))
intplot

intplot.cont<-intresults.cont%>%
  filter(grepl("z.",term))%>%
  mutate(ci.high=estimate+1.64*std.error,
         ci.low=estimate-1.64*std.error)%>%
  ggplot(aes(x=estimate,y=term))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high))+
  geom_vline(xintercept=0,lty="dotted")+
  scale_y_discrete(limits=c("zincrease.climate","zdecrease.climate","zincrease.security","zdecrease.security","zincrease.health","zdecrease.health"),
                   labels = c("zincrease.climate"= "Higher cost,\nclimate change","zdecrease.climate"="Lower cost,\nclimate change", 
                              "zincrease.security"="Higher cost,\nenergy security","zdecrease.security"="Lower cost,\nenergy security",
                              "zincrease.health"="Higher cost,\nhealth","zdecrease.health"="Lower cost,\nhealth")
  ) +
  labs(x="",y="")+
  coord_cartesian(xlim = c(-0.2,0.1)) +
  theme_classic()+
  theme(axis.text.x = element_text(size=11), 
        axis.text.y = element_text(size=11))
intplot.cont

## figure 4: combined results plot #### 
results.combined<-annotate_figure(ggarrange(intplot,mainplot),
                                  left=text_grob("Treatment",size=11,rot=90),
                                  bottom = text_grob("Effect on policy support, relative to control", size = 11))
ggsave(file=paste0(figures,"results_combined.pdf"),height=3.5,width=6.5,results.combined)

## figure S2: combined results plot with ordinal DV
results.combined.cont<-annotate_figure(ggarrange(intplot.cont,mainplot.cont),
                                       left=text_grob("Treatment",size=11,rot=90),
                                       bottom = text_grob("Effect on policy support, relative to control", size = 11))
ggsave(file=paste0(figures,"results_combined_continuousDV.pdf"),height=3.5,width=6.5,
       results.combined.cont)

# minimum detectable effect #### 
## standard error of the average treatment effect for the low-cost treatment: 
## pg. 61 from Gerber and Green 2012 (eqn. 3.6)
## var Y0: 
d2<-d%>%
  filter(z.cost=="no cost info"|z.cost=="high")%>%
  mutate(d.diff=case_when(z.cost=="no cost info"~dv.methanelaw.binary-mean(d$dv.methanelaw.binary[d$z.cost=="no cost info"]),
                          z.cost=="high"~dv.methanelaw.binary-mean(d$dv.methanelaw.binary[d$z.cost=="high"])))
se_y0_hat<-sqrt(sum((d2$d.diff[d2$z.cost=="no cost info"])^2)/(nrow(d2[d2$z.cost=="no cost info",])-1))
se_y1_hat<-sqrt(sum((d2$d.diff[d2$z.cost=="high"])^2)/(nrow(d2[d2$z.cost=="high",])-1))
# sd_y_0_hat <- sd(d$dv.methanelaw.binary[d$z.cost == "no cost info"])
# sd_y_1_hat <- sd(d$dv.methanelaw.binary[d$z.cost == "high"])
# sd_y_2_hat <- sd(d$dv.methanelaw.binary[d$z.cost == "low"])
se_hat_full_study <- sqrt( se_y0_hat^2 / nrow(d2[d2$z.cost=="no cost info",])  + se_y1_hat^2 / nrow(d2[d2$z.cost=="high",]))
mde_hat_full_study <- 2.8 * se_hat_full_study

# results by knowledge ####
table(d$gascomposition.binary)
table(d$z.cost)
## cost regressions by knowledge #### 
d$knowledge<-factor(d$gascomposition.binary,levels=c(0,1),labels=c("low","high"))
w <- weightit(z.cost~block,data=d[d$gascomposition.binary==1,],estimand="ATE",method="ps")
d[d$gascomposition.binary==1,"w.cost.byknowledge"]<-w$weights
w<-weightit(z.cost~block,data=d[d$gascomposition.binary==0,],estimand="ATE",method="ps")
d[d$gascomposition.binary==0,"w.cost.byknowledge"]<-w$weights
h1mod.know1<-tidy(lm_robust(dv.methanelaw.binary~z.cost+country+ideo+gender+age+education_recode+income_recode,
                            weights=w.cost.byknowledge,data=subset(d,gascomposition.binary==1)))%>%
  mutate(knowledge="high")
h1mod.know0<-tidy(lm_robust(dv.methanelaw.binary~z.cost+country+ideo+gender+age+education_recode+income_recode,
                            weights=w.cost.byknowledge,data=subset(d,gascomposition.binary==0)))%>%
  mutate(knowledge="low")
costint<-tidy(lm_robust(dv.methanelaw.binary~z.cost*knowledge+country+ideo+gender+age+education_recode+income_recode,
                          weights=w.cost,data=d))%>%
  filter(grepl(":",term))%>%
  mutate(sig=ifelse(p.value<0.05,1,0))%>%
  select(term,p.value,sig)## neither effect is sig

h1mod.knowledge<-rbind(h1mod.know1,h1mod.know0)

## benefit frames regressions by knowledge #### 
w<-weightit(z.frame~block,estimand="ATE",method="ps",data=d[d$gascomposition.binary==1,])
d[d$gascomposition.binary==1,"w.frame.know"]<-w$weights
w<-weightit(z.frame~block,estimand="ATE",method="ps",data=d[d$gascomposition.binary==0,])
d[d$gascomposition.binary==0,"w.frame.know"]<-w$weights

frameresults.know1<-tidy(lm_robust(dv.methanelaw.binary~z.frame+country+ideo+gender+age+education_recode+income_recode,
                              weights=w.frame.know,data=subset(d,gascomposition.binary==1)))%>%
  mutate(knowledge="high")
frameresults.know0<-tidy(lm_robust(dv.methanelaw.binary~z.frame+country+ideo+gender+age+education_recode+income_recode,
                              weights=w.frame.know,data=subset(d,gascomposition.binary==0)))%>%
  mutate(knowledge="low")
h2mod.knowledge<-rbind(frameresults.know1,frameresults.know0)
frameint<-tidy(lm_robust(dv.methanelaw.binary~z.frame*knowledge+country+ideo+gender+age+education_recode+income_recode,
                         weights=w.frame,data=d))%>%
  mutate(sig=ifelse(p.value<0.05,1,0))%>%
  filter(grepl(":",term))%>%
  select(term,p.value,sig) ## all null 

## interaction effects by knowledge #### 
w<- weightit(z~block,data=d[d$gascomposition.binary==1,],estimand="ATE",method="ps")
d[d$gascomposition.binary==1,"w.full.know"]<-w$weights
w<- weightit(z~block,data=d[d$gascomposition.binary==0,],estimand="ATE",method="ps")
d[d$gascomposition.binary==0,"w.full.know"]<-w$weights

h3amod.know1<-tidy(lm_robust(dv.methanelaw.binary~z+country+ideo+gender+age+education_recode+income_recode,
                  data=subset(d,gascomposition.binary==1),weights=w.full.know))%>%
  mutate(knowledge="high")
h3amod.know0<-tidy(lm_robust(dv.methanelaw.binary~z+country+ideo+gender+age+education_recode+income_recode,
                             data=subset(d,gascomposition.binary==0),weights=w.full.know))%>%
  mutate(knowledge="low")
h3amod.know<-rbind(h3amod.know0,h3amod.know1)

fullint<-tidy(lm_robust(dv.methanelaw.binary~z*knowledge+country+ideo+gender+age+education_recode+income_recode,
                         weights=w.frame,data=d))%>%
  filter(grepl(":",term))%>%
  mutate(sig=ifelse(p.value<0.05,1,0))%>%
  select(term,p.value,sig) ## increase.health is only significant difference

## table S10 of results by knowledge #### 
resultstab.know<-do.call(rbind,list(h3amod.know,h1mod.knowledge,h2mod.knowledge))%>%
  filter(grepl("z",term))
resultstab.know$`BH p value` <- round(p.adjust(resultstab.know$p.value,"BH"),2)
resultstab.know<-resultstab.know%>%
  select(term,estimate,std.error,statistic,p.value,`BH p value`,knowledge)%>%
  mutate(term=case_when(term=="zincrease.security"~"Higher cost, energy security",
                 term=="zdecrease.security"~"Lower cost, energy security",
                 term=="zincrease.climate"~"Higher cost, climate change",
                 term=="zdecrease.climate"~"Lower cost, climate change",
                 term=="zincrease.health"~"Higher cost, health",
                 term=="zdecrease.health"~"Lower cost, health",
                 term=="z.costlow"~ "Lower costs",
                 term=="z.costhigh"~"Higher costs",
                 term=="z.framehealth"~"Health",
                 term=="z.frameclimate"~"Climate change",
                 term=="z.frameenergysecurity"~"Energy security"))
print.xtable(xtable(resultstab.know,caption="{\\bf Experimental effects by knowledge:} The table shows the results from OLS regressions of our treatment indicators for the cost, benefit, and cost x benefit treatments on split samples of respondents defined by their knowledge about methane in the atmosphere. The penultimate column shows p-values after applying the Benjamini-Hochberg correction for multiple hypothesis testing.",
                    label="tab:knowledgeresults"),include.rownames = FALSE,file=paste0(tables,"results_knowledgeresults.tex"))

## Fig. S3: 2-panel effects plot by knowledge #### 
knowplot.right<-rbind(h1mod.knowledge,h2mod.knowledge)%>%
  filter(grepl("z.",term))%>%
  mutate(ci.high=estimate+1.96*std.error,
         ci.low=estimate-1.96*std.error)%>%
  ggplot(aes(x=estimate,y=term,shape=knowledge))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high),position=position_dodge(width=.2))+
  geom_vline(xintercept=0,lty="dotted")+
  theme_classic()+
  scale_y_discrete(
    labels = c("z.costlow" = "Lower costs",
               "z.costhigh" = "Higher costs",
               "z.framehealth" = "Health",
               "z.frameclimate" = "Climate change",
               "z.frameenergysecurity" = "Energy security"))+
  labs(x="",y="")+
  scale_x_continuous(limits=c(-.15,.15))+
  # coord_cartesian(xlim = c(-0.1,0.1)) +
  theme(axis.text.x = element_text(size=11),
        axis.text.y = element_text(size=11))
knowplot.right

knowplot.left<-h3amod.know%>%
  filter(grepl("z.",term))%>%
  mutate(ci.high=estimate+1.96*std.error,
         ci.low=estimate-1.96*std.error)%>%
  ggplot(aes(x=estimate,y=term,shape=knowledge))+
  geom_pointrange(aes(xmin=ci.low,xmax=ci.high),position=position_dodge(width=.2))+
  geom_vline(xintercept=0,lty="dotted")+
  scale_y_discrete(
    limits = c("zincrease.security","zdecrease.security",
               "zincrease.climate","zdecrease.climate",
               "zincrease.health","zdecrease.health"),
    labels = c("Higher cost,\nenergy security","Lower cost,\nenergy security",
               "Higher cost,\nclimate change","Lower cost,\nclimate change", 
               "Higher cost,\nhealth","Lower cost,\nhealth")
  ) +
  labs(x="",y="")+
  scale_x_continuous(limits=c(-.15,.15))+
  # coord_cartesian(xlim = c(-0.1,0.1)) +
  theme(legend.position = "bottom",
        panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        axis.line.x = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.x = element_text(size=11), 
        axis.title.x = element_text(size=13), 
        axis.line.y = element_line(colour = 'black', size=0.5, linetype='solid'),
        axis.text.y = element_text(size=11), 
        axis.title.y = element_text(size=13))
knowplot.left

results.combined.byknowledge<-annotate_figure(ggarrange(knowplot.left,knowplot.right,common.legend=TRUE,
                                                        legend="bottom"),
                                  left=text_grob("Treatment",size=11,rot=90),
                                  bottom = text_grob("Effect on policy support, relative to control", size = 11))
ggsave(file=paste0(figures,"results_combined_byknowledge.pdf"),height=3.5,width=6.5,results.combined.byknowledge)

# exploratory hypotheses #### 
## cost x covariate interactions #### 
eh1amod<-lm_robust(dv.methanelaw.binary~z.cost*country+ideo+gender+age+education_recode+income_recode,
                   weights=w.cost,data=d)

eh1bmod<-lm_robust(dv.methanelaw.binary~z.cost*ideo+country+gender+age+education_recode+income_recode,
                   weights=w.cost,data=d)

eh1cmod<-lm_robust(dv.methanelaw.binary~z.cost*age+country+ideo+gender+education_recode+income_recode,
                   weights=w.cost,data=d)

eh1dmod<-lm_robust(dv.methanelaw.binary~z.cost*education_recode+ideo+gender+age,#+income+education,
                   weights=w.cost,data=d)

eh1emod<-lm_robust(dv.methanelaw.binary~z.cost*income_recode+ideo+gender+age,#+income+education,
                   weights=w.cost,data=d)

## table S11: conditional effects of cost information #### 
texreg(list(eh1amod,eh1bmod, eh1cmod,eh1dmod,eh1emod),
       custom.model.names = c("x countries","x ideologies","x ages",
                              "x education","x incomes"),
       include.ci=FALSE,caption="Conditional effects of cost information",
       custom.coef.names=c("(Intercept)","High cost","Low cost","Germany","Italy","Poland",
                           "Left","Right","Male","Age",
                           "Secondary ed or less","Income",
                           "High cost x Germany","Low cost x Germany",
                           "High cost x Italy","Low cost x Italy",
                           "High cost x Poland","Low cost x Poland",
                           "High cost x Left","Low cost x Left",
                           "High cost x Right","Low cost x Right",
                           "High cost x Age","Low cost x Age",
                           "High cost x Ed","Low cost x Ed",
                           "High cost x Income","Low cost x Income"),
       fontsize="footnotesize",single.row=TRUE,
       file=paste0(tables,"results_costinteractions.tex"),
       label="tab:costinteractions")

## benefit x covariate interactions #### 
eh2amod<-lm_robust(dv.methanelaw.binary~z.frame*country+ideo+gender+age+education_recode+income_recode,
                   weights=w.frame,data=d)

eh2bmod<-lm_robust(dv.methanelaw.binary~z.frame*ideo+country+gender+age+education_recode+income_recode,
                   weights=w.frame,data=d)

eh2cmod<-lm_robust(dv.methanelaw.binary~z.frame*age+country+ideo+gender+education_recode+income_recode,
                   weights=w.frame,data=d)

eh2dmod<-lm_robust(dv.methanelaw.binary~z.frame*education_recode+country+ideo+gender+age+income_recode,weights=w.frame,data=d) ##  add this and next regression in when we have education and income data coded.

eh2emod<-lm_robust(dv.methanelaw.binary~z.frame*income_recode+country+ideo+gender+education_recode+age,weights=w.frame,data=d)


## table S12: conditional effects of benefit framing #### 
texreg(list(eh2amod,eh2bmod,eh2cmod,eh2dmod,eh2emod
),
custom.model.names = c("x countries","x ideologies","x ages","x education","x incomes"
),
custom.coef.names=c("(Intercept)","Climate frame","Security frame","Health frame","Germany","Italy","Poland",
                    "Left","Right","Male","Age",
                    "Secondary ed or less","Income",
                    "Climate frame x Germany","Security frame x Germany",
                    "Health frame x Germany",
                    "Climate frame x Italy","Security frame x Italy",
                    "Health frame x Italy",
                    "Climate frame x Poland","Security frame x Poland",
                    "Health frame x Poland",
                    "Climate frame x Left","Security frame x Left",
                    "Health frame x Left",
                    "Climate frame x Right","Security frame x Right",
                    "Health frame x Left",
                    "Climate frame x Age","Security frame x Age",
                    "Health frame x Age",
                    "Climate frame x Ed","Security frame x Ed",
                    "Health frame x Ed",
                    "Climate frame x Income","Security frame x Income",
                    "Health frame x"),
include.ci=FALSE,caption="Conditional effects of benefit framing",
fontsize="footnotesize",single.row=TRUE,
file=paste0(tables,"results_frameinteractions.tex"),
label="tab:frameinteractions",float.pos="!htpb")

## conditional effects of cost-sensitivity of responses to benefit frames #### 
eh3amod<-lm_robust(dv.methanelaw.binary~z*country+ideo+gender+age+education_recode+income_recode,
                   weights=w.full,data=d)

eh3bmod<-lm_robust(dv.methanelaw.binary~z*ideo+country+gender+age+education_recode+income_recode,
                   weights=w.full,data=d)

eh3cmod<-lm_robust(dv.methanelaw.binary~z*age+country+ideo+gender+education_recode+income_recode,
                   weights=w.full,data=d)

eh3dmod<-lm_robust(dv.methanelaw.binary~z*education_recode+country+ideo+gender+age+income_recode,weights=w.full,data=d) 

eh3emod<-lm_robust(dv.methanelaw.binary~z*income_recode+country+ideo+gender+education_recode+age,weights=w.full,data=d)

## table S13: conditional effects of cost sensitivity of benefit frames #### 
texreg(list(eh3amod,eh3bmod,eh3cmod,eh3dmod,eh3emod
),
custom.model.names = c("x countries","x ideologies","x ages","x education","x incomes"
),caption="Conditional effects of cost sensitivity of benefit frames",
float.pos="!htpb",
include.ci=FALSE,fontsize="tiny",single.row=TRUE,
custom.coef.names=c("Intercept","Low cost x climate","Low cost x health","Low cost x security",
                    "High cost x climate","High cost x health","High cost x security",
                    "Germany","Italy","Poland","Left","Right","Male","Age","Secondary ed or less","Income",
                    "Low cost x climate x Germany","Low cost x health x Germany","Low cost x security x Germany",
                    "High cost x climate x Germany","High cost x health x Germany","High cost x security x Germany",
                    "Low cost x climate x Italy","Low cost x health x Italy","Low cost x security x Italy",
                    "High cost x climate x Italy","High cost x health x Italy","High cost x security x Italy",
                    "Low cost x climate x Poland","Low cost x health x Poland","Low cost x security x Poland",
                    "High cost x climate x Poland","High cost x health x Poland","High cost x security x Poland",
                    "Low cost x climate x Left","Low cost x health x Left","Low cost x security x Left",
                    "High cost x climate x Left","High cost x health x Left","High cost x security x Left",
                    "Low cost x climate x Right","Low cost x health x Right","Low cost x security x Right",
                    "High cost x climate x Right","High cost x health x Right","High cost x security x Right",
                    "Low cost x climate x age","Low cost x health x age","Low cost x security x age",
                    "High cost x climate x age","High cost x health x age","High cost x security x age",
                    "Low cost x climate x ed","Low cost x health x ed","Low cost x security x ed",
                    "High cost x climate x ed","High cost x health x ed","High cost x security x ed",
                    "Low cost x climate x income","Low cost x health x income","Low cost x security x income",
                    "High cost x climate x income","High cost x health x income","High cost x security x income"),
file=paste0(tables,"results_fullinteractions.tex"),label="tab:fullinteractions")

